MODULE init
  
  USE definitions
  
  IMPLICIT NONE
  
  
CONTAINS
  
  SUBROUTINE init_relpos
    
    USE definitions
    
    IMPLICIT NONE
    
    relpos(1,1) = -1
    relpos(1,2) = +1
    relpos(1,3) = 0
    relpos(1,4) = 0
    relpos(1,5) = -1
    relpos(1,6) = +1
    relpos(1,7) = -1
    relpos(1,8) = +1
    relpos(2,1) = 0
    relpos(2,2) = 0
    relpos(2,3) = -1
    relpos(2,4) = +1
    relpos(2,5) = -1
    relpos(2,6) = -1
    relpos(2,7) = +1
    relpos(2,8) = +1
    
  END SUBROUTINE init_relpos
  
  SUBROUTINE init_grid
    
    
    USE definitions
    
    IMPLICIT NONE
    
    INTEGER :: i, j, x, y		!Counters for Do loops
    REAL :: Slabnum
    
    IF(load == 0) THEN
       WRITE(*,*) 'initialize Sandpit to', xsize, 'times', ysize
       
       ALLOCATE(Sandpit(1:xsize,1:ysize))
       
       DO i=1, xsize
          DO j=1, ysize
             
             Sandpit(i,j)%Slabnumber=initialhight
             
          END DO
       END DO
       
    ELSEIF(load == 1) THEN
       
       WRITE(*,*) 'load state from final.dat'
       WRITE(*,*) 'parameters are used from ast.in'
       OPEN(unit = in1, file = 'final.dat',status = 'old', action='read')
       
       READ(in1, *) xsize, ysize
       WRITE(*,*) 'width of array is', xsize, 'times', ysize
       
       ALLOCATE(Sandpit(1:xsize,1:ysize))
       
       DO i=1,xsize
          DO j=1,ysize
             
             READ(in1, *) x, y, Slabnum
             
             Sandpit(x,y)%Slabnumber = Slabnum
             
          ENDDO
       ENDDO
       CLOSE(unit = in1)
    ELSEIF(load == 2) THEN
       ALLOCATE(Sandpit(1:xsize,1:ysize))
       DO i = 1,xsize
          DO j = 1,ysize
             IF(i<10) THEN
                Sandpit(i,j)%slabnumber = initialhight
             ELSEIF(i>=10) THEN
                Sandpit(i,j)%slabnumber = 1
             ENDIF
          ENDDO
       ENDDO
    ENDIF

    ALLOCATE(position(1:xsize*ysize))
    
    WRITE(*,*) xsize*ysize

    DO i=1,xsize
       DO j=1,ysize
          position((i-1)*ysize+j)%entry(1) = i
          position((i-1)*ysize+j)%entry(2) = j
       ENDDO
    ENDDO

  END SUBROUTINE init_grid
  
  SUBROUTINE init_rand(seed)
    
    IMPLICIT NONE
    
    INTEGER, INTENT(inout) :: seed
    
    WRITE(*,*) 'Seed for random numbers generator is', seed
    
    CALL random_seed(seed)
    
    WRITE(*,*) 'Output of random generator initialisation is', seed
    
  END SUBROUTINE init_rand
  
  SUBROUTINE init_counter
    
    USE definitions
    
    IMPLICIT NONE

    CALL SYSTEM_CLOCK(count1, countrate, countmax)
    
    avalcount = 0
    eshadowcount = 0
    dshadowcount = 0
    empty = 0
    avhightbuf = 0
    
  END SUBROUTINE init_counter
  
  SUBROUTINE init_output
    
    USE definitions

    IMPLICIT NONE

    CALL SYSTEM('rm -r data')
    CALL SYSTEM('mkdir data')

    ALLOCATE(max_hight(xsize))
    max_hight = 0

    ALLOCATE(local_max_hight(xsize,ysize))
    local_max_hight = 0

    OPEN(unit = out4, file = 'dynamics.dat')
    IF(slices == 1) CALL init_slice_diagnostics
    IF(full == 1) CALL init_full_diagnostics

    CALL SYSTEM('sh initdata.sh')
    
  END SUBROUTINE init_output

  SUBROUTINE final_output

    USE definitions
    USE inout
    
    IMPLICIT NONE

    INTEGER :: i, j, k, l, m, covered_cells
    INTEGER, DIMENSION(2) :: pos
    REAL    :: minhight, maxhight
    
    OPEN(unit = out21, file='diag.dat', access='APPEND')
    
    WRITE(*,*) 'calculating diagnostic values'
    
    IF(source == 1) THEN
       maxhight = sourceflow
       minhight = sourceflow
    ELSEIF(initialhight>0 .AND. source == 0) THEN
       maxhight = initialhight
       minhight = initialhight
    ELSEIF(initialhight>0 .AND. source == 1) THEN
       maxhight = initialhight + sourceflow
       minhight = initialhight + sourceflow
    ENDIF

    covered_cells = 0

    DO k=2,ysize-1
       DO j=2,xsize-1
          pos(1) = j
          pos(2) = k
          IF(sandpit(j,k)%slabnumber>maxhight) maxhight = sandpit(j,k)%slabnumber
          IF(sandpit(j,k)%slabnumber<minhight) minhight = sandpit(j,k)%slabnumber
          IF(sandpit(j,k)%slabnumber>0) covered_cells = covered_cells + 1
       ENDDO
    ENDDO
    WRITE(*,*)'writing diagnostics to diag.dat'
    CALL SYSTEM_CLOCK(count2, countrate, countmax)
    WRITE(*,*) REAL((count2-count1))/REAL(countrate), sourceflow, minhight, maxhight, REAL(covered_cells)/REAL(xsize*ysize)
    WRITE(out21, '(4I8, A, 9E16.8)') xsize, ysize, steps, &
         INT(REAL((count2-count1))/REAL(countrate)), char(9), &
         sourceflow, L0, E0, minhight, maxhight, REAL(covered_cells)/REAL(xsize*ysize), dwflux,&
         daflux, initialhight
         
    
    CLOSE(out21)

    CALL write_envelope
    CALL write_final_state			!Write final slab hight to final.dat
    CALL write_parameters			!Write simulation Parameters to File
    CALL fft
    IF(slices == 1) CALL finalize_slice_diagnostics
    IF(full == 1) CALL finalize_full_diagnostics
    CALL SYSTEM('sh plotdiag.sh')
    CALL SYSTEM('sh movedata.sh')		!Move output to data folder
    CALL SYSTEM('sh slabhight.sh')
    
    CLOSE(unit = out4)
    
    DEALLOCATE(Sandpit)
    DEALLOCATE(position)
    DEALLOCATE(max_hight)
    DEALLOCATE(local_max_hight)

  END SUBROUTINE final_output

  !------------------------------------------------------------
  !Initialize slice diags: open files, initialize counters and other stuff...
  !------------------------------------------------------------
  SUBROUTINE init_slice_diagnostics

    USE definitions
    
    IMPLICIT NONE

    INTEGER :: irec

    ALLOCATE(sliced_transport_length(0:steps+1,1:ysize))
    ALLOCATE(sliced_erosion_hight(0:steps+1,1:ysize))
    ALLOCATE(sliced_average_hight(0:steps+1,1:ysize))
    ALLOCATE(sliced_average_flux(0:steps+1,1:ysize))
    ALLOCATE(sliced_covered_cells(0:steps+1,1:ysize))

  END SUBROUTINE init_slice_diagnostics
  !------------------------------------------------------------

  !------------------------------------------------------------
  !finalize slice diagnostics
  !------------------------------------------------------------
  SUBROUTINE finalize_slice_diagnostics

    USE definitions

    IMPLICIT NONE

    INTEGER :: irec, t, y
    REAL, DIMENSION(ysize) :: covered_cells
    
    
       INQUIRE(iolength = irec) sliced_transport_length
       OPEN(unit=out6, file='sliced_transport_length.dat', form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out6, rec = 1) sliced_transport_length

       INQUIRE(iolength = irec) sliced_erosion_hight
       OPEN(unit=out7, file='sliced_erosion_hight.dat', form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out7, rec = 1) sliced_erosion_hight

       INQUIRE(iolength = irec) sliced_average_hight
       OPEN(unit=out8, file='sliced_average_hight.dat', form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out8, rec = 1) sliced_average_hight

       INQUIRE(iolength = irec) sliced_average_flux
       OPEN(unit=out9, file='sliced_average_flux.dat', form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out9, rec = 1) sliced_average_flux

       INQUIRE(iolength = irec) sliced_covered_cells
       OPEN(unit=out10, file='sliced_covered_cells.dat', form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out10, rec = 1) sliced_covered_cells
       CLOSE(unit = out6)
       CLOSE(unit = out7)
       CLOSE(unit = out8)
       CLOSE(unit = out9)
       CLOSE(unit = out10)

       covered_cells = 0
       
       DO y=1,ysize
          DO t=tconv,steps 
             covered_cells(y) = covered_cells(y) + sliced_covered_cells(t,y)
          ENDDO
       ENDDO
       
       covered_cells = covered_cells/REAL(steps-tconv)

       OPEN(unit = out11, file='mean_covered_cells.dat')
WRITE(*,*) covered_cells
       DO y=1,ysize
          WRITE(out11, '( I4, E16.8 )') y, covered_cells(y)
       ENDDO

       CLOSE(out11)

    DEALLOCATE(sliced_transport_length)
    DEALLOCATE(sliced_erosion_hight)
    DEALLOCATE(sliced_average_hight)
    DEALLOCATE(sliced_average_flux)
    DEALLOCATE(sliced_covered_cells)

  END SUBROUTINE finalize_slice_diagnostics
  !------------------------------------------------------------

  SUBROUTINE init_full_diagnostics

    USE definitions

    IMPLICIT NONE

    mintl=0
    maxtl=0
    mineh=0
    maxeh=0
    mindh=0
    maxdh=0
    minaf=0
    maxaf=0
    minvf=0
    maxvf=0
    minsb=0
    maxsb=0

    ALLOCATE(local_transport_length(1:xsize,1:ysize))
    ALLOCATE(local_eroded_hight(1:xsize,1:ysize))
    ALLOCATE(local_deposited_hight(1:xsize,1:ysize))
    ALLOCATE(local_aeolian_flux(1:xsize,1:ysize))
    ALLOCATE(local_avalanche_flux(1:xsize,1:ysize))
    ALLOCATE(local_sediment_balance(1:xsize,1:ysize))

  END SUBROUTINE init_full_diagnostics

  SUBROUTINE finalize_full_diagnostics

    USE definitions

    IMPLICIT NONE
    

    OPEN(unit=out4, file = 'parameters.dat', status = 'old', position = 'append')

    WRITE(out4, '(A, E16.8, A, A5)') '#', mintl, char(9), 'mintl'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxtl, char(9), 'maxtl'
    WRITE(out4, '(A, E16.8, A, A5)') '#', mineh, char(9), 'mineh'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxeh, char(9), 'maxeh'
    WRITE(out4, '(A, E16.8, A, A5)') '#', mindh, char(9), 'mindh'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxdh, char(9), 'maxdh'
    WRITE(out4, '(A, E16.8, A, A5)') '#', minaf, char(9), 'minaf'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxaf, char(9), 'maxaf'
    WRITE(out4, '(A, E16.8, A, A5)') '#', minvf, char(9), 'minvf'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxvf, char(9), 'maxvf'
    WRITE(out4, '(A, E16.8, A, A5)') '#', minsb, char(9), 'minsb'
    WRITE(out4, '(A, E16.8, A, A5)') '#', maxsb, char(9), 'maxsb'

    CLOSE(unit = out4)

    CALL SYSTEM('sh plotfulldiag.sh')

    DEALLOCATE(local_transport_length)
    DEALLOCATE(local_eroded_hight)
    DEALLOCATE(local_deposited_hight)
    DEALLOCATE(local_aeolian_flux)
    DEALLOCATE(local_avalanche_flux)
    DEALLOCATE(local_sediment_balance)

    END SUBROUTINE finalize_full_diagnostics
  
END MODULE init


